#####################################################################
#Program for fit our propose model considering \alpha_1 and \alpha_2 equal to a constant (Method 4)
#####################################################################

misclass_a1_a2=function (form,y1=y1,y2=y2, a1 = 0, a2 = 0, bmat = 0, print.summary = TRUE, 
                         data = NULL) # form= y~x
{
  probit <- glm(form, family = binomial(link = "probit"), data = data)
  cat("Standard Probit Estimates", "\n")
  print(summary(probit))
  xmat <- model.matrix(form, data = data)
  if (identical(bmat, 0)) {
    bmat <- probit$coef
  }
  nk = length(bmat)
  y <- model.frame(form, data = data)[, 1]
  a1 = ifelse(identical(a1, 0), 0.001, a1)
  a2 = ifelse(identical(a2, 0), 0.001, a2)
  bstart <- c(bmat, qnorm(a1), qnorm(a2))
  logl <- function(b) {
    y00 <- (1 - y1)*(1 - y2)
    y01 <- (1 - y1)*y2
    y10 <- y1*(1 - y2) 
    y11<- y1*y2
    
    p00<- 1- pnorm(xmat %*% b[1:nk])+ pnorm(b[nk + 1])* pnorm(b[nk + 2])* pnorm(xmat %*% b[1:nk])
    logprob00 <- log(p00)
    p01<- pnorm(b[nk + 1])*(1- pnorm(b[nk + 2]))* pnorm(xmat %*% b[1:nk])
    logprob01 <- log(p01)
    p10<- (1-pnorm(b[nk + 1]))* pnorm(b[nk + 2])* pnorm(xmat %*% b[1:nk])
    logprob10 <- log(p10)
    p11<- (1-pnorm(b[nk + 1]))*(1- pnorm(b[nk + 2]))* pnorm(xmat %*% b[1:nk])
    logprob11 <- log(p11)
    y00t <- t(y00)
    y01t <- t(y01)
    y10t <- t(y10)
    y11t <- t(y11)
    out <- - sum(y00t%*%logprob00 + y01t%*%logprob01 + y10t%*%logprob10 +
                   y11t%*%logprob11)           
    return(out)
  }
  nlfit<- optim( bstart,logl,method="L-BFGS-B",control=list(maxit=500), hessian = TRUE)
  vmat <- solve(nlfit$hessian)
  stderr <- sqrt(abs(diag(vmat)))
  tval <- nlfit$par/stderr
  pval <- 2 * (1 - pnorm(abs(tval)))
  nall <- length(nlfit$par)
  nk1 = nk + 1
  amat <- pnorm(nlfit$par[nk1:nall])
  tmat <- tval[nk1:nall]
  smat <- abs(amat/tmat)
  nvect <- c("a1", "a2")
  names(nlfit$par[nk1:nall]) <- nvect
  names(amat) <- nvect
  names(tmat) <- nvect
  names(smat) <- nvect
  cname <- c("Estimate", "Std. Error", "z-value", "Pr(>|z|)")
  if (print.summary == TRUE) {
    outmat <- cbind(nlfit$par, stderr, tval, pval)
    rownames(outmat) <- c(names(probit$coef), paste("qnorm(", 
                                                    nvect, ")", sep = ""))
    colnames(outmat) <- cname
    cat("nlm estimates", "\n")
    cat(" ", "\n")
    print(round(outmat, 5))
    cat(" ", "\n")
    cat("Estimated misclassification probabilities:", "\n")
    outmat <- cbind(amat, smat, tmat, 2 * (1 - pnorm(abs(tmat))))
    rownames(outmat) <- nvect
    colnames(outmat) <- cname
    print(round(outmat, 5))
    cat(" ", "\n")
  }
  
  out1 <- list(nlfit$par, stderr, vmat, nlfit$convergence, 
               -nlfit$value, amat[1], amat[2],smat)
  names(out1) <- c( "estimate", "stderr", "vmat","convergence", "minimum","a1", "a2", "smat")
  return(out1)  
}
